C======================================================================      
C   This subroutine sets up the r-mesh by Simpson rule.
C======================================================================      
      subroutine  gridsm (qmax, rmax, gridr, w, nr, nload)
      dimension  gridr(nload), w(nload) 
C
C  INPUT:
C  -----
C   qmax - the biggest momentum what could be calculated by this mesh (a.u.).
C   rmax - the largest "r" in the mesh.
C   nload- the load dimention of arries gridr and  w
C
C  OUTPUT:
C  ------
C   nr - total number of points in the mesh.
C   gridr(i) - knots.   
C   w(i)     - weigths. 
C
C--------------------------------------------------------------------------
C
C     Do some check 
C
      step = 3.14 / 6.0 / qmax
      if (step .le. 0.0) then
         print*, 'Step is wrong', step
         stop    'Step is wrong'
      end if 
      npoints = nint(rmax / step)
      if (npoints .gt. nload) then
         print*, 'Increase number of points (nload)', npoints
         stop    'Increase number of points'
      end if 
C
C     Set up r-mesh 
C
      dr = step
      r = 0.0
      do i = 1, npoints
         r = r + dr
         gridr(i) = r 
         w(i) = float(mod(i, 2) * 2 + 2) * dr / 3.0
      end do 
      nr = npoints
      return
      end




